///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file agp.cpp
///
/// \brief   contains implementations of functions related to the antisymmetrized geminal power
///          wavefunction
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include <src/pcps.h>
#include <src/agp.h>
#include <src/matrix_algebra.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   checks that the wavefunction's data arrays are of the expected size for working with
///          AGP and pairing matrix variables
///
/// \param[in]     wfn      the wavefunction
/// \param[in]     nap      number of alpha particles
/// \param[in]     nas      number of alpha lattice sites
/// \param[in]     use_pm   whether to use the pairing matrix or AGP weights and orbital coefficients
///
/// \return true
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class SCALAR> bool pcps::agp::check_sizes(const pcps::Wavefunction<SCALAR> & wfn,
                                                    const int nap,
                                                    const int nas,
                                                    const bool use_pm)
{
  assert(  use_pm || wfn.agp_weights().size() == nas );
  assert(  use_pm || wfn.restricted_orb_coeffs().size() == nas*nas );
  assert( !use_pm || wfn.pairing_matrix().size() == nas*nas );
  return true;
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   populates the arrays of occupied alpha and beta sites for the given configuration
///
/// \param[in]     nap       number of alpha particles
/// \param[in]     nas       number of alpha lattice sites
/// \param[in]     config    lattice configuration, size 2*nas
/// \param[out]    a_sites   work space to store the occupied alpha sites, size nap
/// \param[out]    b_sites   work space to store the occupied beta  sites, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void pcps::agp::populate_site_lists(const int nap,
                                    const int nas,
                                    const int * const config,
                                    int * const a_sites,
                                    int * const b_sites)
{
  int a_count = 0;
  int b_count = 0;
  for (int i = 0; i < 2*nas; i+=2) {
    if (config[i] == 1)
      a_sites[a_count++] = i/2;
    if (config[i+1] == 1)
      b_sites[b_count++] = i/2;
  }
  assert(a_count == nap);
  assert(b_count == nap);
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   builds the occupied pairing matrix
///
/// \param[in]     nap       number of alpha particles
/// \param[in]     nas       number of alpha lattice sites
/// \param[in]     gamma     the all-sites pairing matrix, size nas*nas
/// \param[in]     a_sites   the occupied alpha sites, size nap
/// \param[in]     b_sites   the occupied beta  sites, size nap
/// \param[out]    mat       on exit, the occupied pairing matrix, size nap*nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class SCALAR> void pcps::agp::build_occupied_matrix(const int nap,
                                                              const int nas,
                                                              const SCALAR * const gamma,
                                                              const int * const a_sites,
                                                              const int * const b_sites,
                                                              SCALAR * const mat)
{
  for (int i = 0; i < nap; i++)
  for (int j = 0; j < nap; j++)
    mat[i*nap+j] = gamma[ a_sites[i] * nas + b_sites[j] ];
}

//-------------------------------------------------------------------------------
// pcps::agp::init_matrices -- Initializes the all-sites pairing matrix and the
//                             inverse of the occupied pairing matrix.
//                             Returns the value of the agp matrix's determinant.
//-------------------------------------------------------------------------------

                                                                       // exit status  array size    description
template <class SCALAR>
void pcps::agp::init_matrices(const pcps::Wavefunction<SCALAR> & wfn,  // unchanged                  the wavefunction object
                              const int nap,                           // unchanged                  number of alpha particles
                              const int nas,                           // unchanged                  number of alpha lattice sites
                              const int * const config,                // unchanged      2*nas       lattice configuration
                              SCALAR & det_value,                      // overwritten                on exit, the value of the agp matrix's determinant
                              int * const a_sites,                     // overwritten    nap         on exit, the occupied alpha sites (numbered 1 to nas)
                              int * const b_sites,                     // overwritten    nap         on exit, the occupied beta sites (numbered 1 to nas)
                              SCALAR * const gamma,                    // overwritten    nas*nas     on exit, the gamma matrix
                              SCALAR * const inv_mat,                  // overwritten    nap*nap     on exit, the inverse of the agp matrix
                              SCALAR * const agp_mat,                  // overwritten    nap*nap     work space for computing the agp matrix
                              int * const ipiv,                        // overwritten    nap         work space for pivoting during the LU decomposition
                              int * const inv_perm,                    // overwritten    nap         work space for inverting the LU pivoting
                              const bool use_pm)              // unchanged                  whether to use the wavefunction's pairing matrix
{

  // ensure the wavefunction's data vectors are the expected size
  assert( pcps::agp::check_sizes(wfn, nap, nas, use_pm) );

  // compute the all-sites pairing matrix
  if (use_pm) {
    pcps::xcopy(nas*nas, &wfn.pairing_matrix()[0], 1, gamma, 1);
  } else {
    for (int p = 0; p < nas; p++)
    for (int q = 0; q < nas; q++) {
      SCALAR & x = gamma[p*nas+q];
      x = pcps::zero<SCALAR>();
      for (int i = 0; i < nas; i++)
        x += wfn.agp_weights()[i] * wfn.restricted_orb_coeffs()[i*nas+p] * wfn.restricted_orb_coeffs()[i*nas+q];
    }
  }

  // populate lists of occupied alpha and beta sites
  pcps::agp::populate_site_lists(nap, nas, config, a_sites, b_sites);

  // build the occupied pairing matrix whose inverse and determinant we need
  pcps::agp::build_occupied_matrix(nap, nas, gamma, a_sites, b_sites, inv_mat);

  // compute the inverse matrix and determinant using an LU decomposition
  pcps::matrix_inverse_lu(nap, det_value, inv_mat, agp_mat, ipiv, inv_perm);

}

// explicitly instantiate the function
template void pcps::agp::init_matrices(const pcps::Wavefunction<double> & wfn,
                                       const int nap,
                                       const int nas,
                                       const int * const config,
                                       double & det_value,
                                       int * const a_sites,
                                       int * const b_sites,
                                       double * const gamma,
                                       double * const inv_mat,
                                       double * const agp_mat,
                                       int * const ipiv,
                                       int * const inv_perm,
                                       const bool use_pm);
template void pcps::agp::init_matrices(const pcps::Wavefunction<std::complex<double> > & wfn,
                                       const int nap,
                                       const int nas,
                                       const int * const config,
                                       std::complex<double> & det_value,
                                       int * const a_sites,
                                       int * const b_sites,
                                       std::complex<double> * const gamma,
                                       std::complex<double> * const inv_mat,
                                       std::complex<double> * const agp_mat,
                                       int * const ipiv,
                                       int * const inv_perm,
                                       const bool use_pm);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the value of the AGP determinant for the supplied configuration using an
///          LU decomposition
///
/// \param[in]     wfn       the wavefunction
/// \param[in]     nap       number of alpha particles
/// \param[in]     nas       number of alpha lattice sites
/// \param[in]     config    lattice configuration, size 2*nas
/// \param[in]     gamma     the all-sites pairing matrix, size nas*nas
/// \param[out]    a_sites   work space to store the occupied alpha sites, size nap
/// \param[out]    b_sites   work space to store the occupied beta  sites, size nap
/// \param[out]    work_mat  work space to store the occupied agp matrix, size nap*nap
/// \param[out]    ipiv      work space for pivoting during the LU decomposition, size nap
/// \param[in]     use_pm    whether to use the pairing matrix or AGP weights and orbital coefficients
///
/// \return the determinant of the occupied AGP matrix in mantissa/exponent format
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class SCALAR> pcps::MantExp<SCALAR> pcps::agp::get_det_value(const pcps::Wavefunction<SCALAR> & wfn,
                                                                       const int nap,
                                                                       const int nas,
                                                                       const int * const config,
                                                                       const SCALAR * const gamma,
                                                                       int * const a_sites,
                                                                       int * const b_sites,
                                                                       SCALAR * const work_mat,
                                                                       int * const ipiv,
                                                                       const bool use_pm)
{

  // ensure the wavefunction's data vectors are the expected size
  assert( pcps::agp::check_sizes(wfn, nap, nas, use_pm) );

  // populate lists of occupied alpha and beta sites
  pcps::agp::populate_site_lists(nap, nas, config, a_sites, b_sites);

  // build the occupied pairing matrix whose determinant we desire
  pcps::agp::build_occupied_matrix(nap, nas, gamma, a_sites, b_sites, work_mat);

  // return the matrix's determinant
  return pcps::matrix_determinant_lu(nap, work_mat, ipiv);

}

// explicitly instantiate the function
template pcps::MantExp<double> pcps::agp::get_det_value(const pcps::Wavefunction<double> & wfn,
                                                        const int nap,
                                                        const int nas,
                                                        const int * const config,
                                                        const double * const gamma,
                                                        int * const a_sites,
                                                        int * const b_sites,
                                                        double * const work_mat,
                                                        int * const ipiv,
                                                        const bool use_pm);
template pcps::MantExp<std::complex<double> > pcps::agp::get_det_value(const pcps::Wavefunction<std::complex<double> > & wfn,
                                                                       const int nap,
                                                                       const int nas,
                                                                       const int * const config,
                                                                       const std::complex<double> * const gamma,
                                                                       int * const a_sites,
                                                                       int * const b_sites,
                                                                       std::complex<double> * const work_mat,
                                                                       int * const ipiv,
                                                                       const bool use_pm);

//-------------------------------------------------------------------------------
// pcps::agp::det_ratio -- Computes the ratio of the AGP determinants for two
//                         configurations that differ at the specified sites.
//                         Returns   det(new config) / det(old config)
//-------------------------------------------------------------------------------

                                                            // exit status  array size    description
template <class SCALAR>
SCALAR pcps::agp::det_ratio(const int nap,                  // unchanged                  number of alpha particles
                            const int nas,                  // unchanged                  number of alpha lattice sites
                            const int nds,                  // unchanged                  the number of sites that differ between the configurations
                            const int * const diff_sites,   // unchanged      nds         the indices of the differing sites ( 0 <= index < 2*nas )
                            const int * const new_config,   // unchanged      2*nas       new lattice configuration
                            const int * const old_config,   // unchanged      2*nas       old lattice configuration
                            const SCALAR * const gamma,     // unchanged      nas*nas     the gamma matrix
                            const SCALAR * const inv_mat,   // unchanged      nap*nap     the inverse of the agp matrix
                            const int * const a_sites,      // unchanged      nap         the occupied alpha sites ( numbered 0 to nas-1 )
                            const int * const b_sites,      // unchanged      nap         the occupied beta sites  ( numbered 0 to nas-1 )
                            SCALAR * const wv1,             // overwritten    nap         work space
                            SCALAR * const wv2,             // overwritten    nap         work space
                            SCALAR * const wv3)             // overwritten    nap         work space
{

  // check that the number of sites affected by the move is supported
  assert(nds == 0 || nds == 2 || nds == 4);

  // return one if no sites were changed
  if (nds == 0)
    return pcps::unity<SCALAR>();

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int n_a_diff, n_b_diff;
  int a_new_occ_sites[2];
  int b_new_occ_sites[2];
  int a_old_occ_sites[2];
  int b_old_occ_sites[2];
  pcps::agp::get_diff_occs(nds, diff_sites, new_config, old_config, n_a_diff, n_b_diff,
                           a_new_occ_sites, b_new_occ_sites, a_old_occ_sites, b_old_occ_sites);

  // create the return value
  SCALAR retval;

  // single alpha move
  if ( n_a_diff == 2 && n_b_diff == 0 ) {

    // get the difference between the new and the old rows
    pcps::agp::row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, wv1);

    // determine which column of the inverse matrix we need
    const int col = pcps::agp::find_index(nap, a_old_occ_sites[0], a_sites);

    // compute the ratio of the new and old determinants by multiplying the row difference by the corresponding column of the inverse matrix
    retval = pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + col, nap, wv1, 1);

  // single beta move
  } else if ( n_a_diff == 0 && n_b_diff == 2 ) {

    // get the difference between the new and the old columns
    pcps::agp::col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, wv1);

    // determine which row of the inverse matrix we need
    const int row = pcps::agp::find_index(nap, b_old_occ_sites[0], b_sites);

    // compute the ratio of the new and old determinants by multiplying the column difference by the corresponding row of the inverse matrix
    retval = pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + nap*row, 1, wv1, 1);

    // double alpha move
  } else if ( n_a_diff == 4 && n_b_diff == 0 ) {

    // get the differences between the new and the old rows
    pcps::agp::row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, wv1);
    pcps::agp::row_diff(nap, nas, a_new_occ_sites[1], a_old_occ_sites[1], gamma, b_sites, wv2);

    // determine which columns of the inverse matrix we need
    const int col1 = pcps::agp::find_index(nap, a_old_occ_sites[0], a_sites);
    const int col2 = pcps::agp::find_index(nap, a_old_occ_sites[1], a_sites);
    assert( col1 != col2 );

    // compute the ratio of the new and old determinants
    retval =   ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + col1, nap, wv1, 1) )
              * ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + col2, nap, wv2, 1) )
              - pcps::xdot(nap, inv_mat + col1, nap, wv2, 1) * pcps::xdot(nap, inv_mat + col2, nap, wv1, 1);

  // double beta move
  } else if ( n_a_diff == 0 && n_b_diff == 4 ) {

    // get the differences between the new and the old columns
    pcps::agp::col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, wv1);
    pcps::agp::col_diff(nap, nas, b_new_occ_sites[1], b_old_occ_sites[1], gamma, a_sites, wv2);

    // determine which row of the inverse matrix we need
    const int row1 = pcps::agp::find_index(nap, b_old_occ_sites[0], b_sites);
    const int row2 = pcps::agp::find_index(nap, b_old_occ_sites[1], b_sites);
    assert( row1 != row2 );

    // compute the ratio of the new and old determinants
    retval =   ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + nap*row1, 1, wv1, 1) )
              * ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + nap*row2, 1, wv2, 1) )
              - pcps::xdot(nap, inv_mat + nap*row1, 1, wv2, 1) * pcps::xdot(nap, inv_mat + nap*row2, 1, wv1, 1);

  // alpha beta move
  } else if ( n_a_diff == 2 && n_b_diff == 2 ) {

    // get the differences between the new and the old rows
    pcps::agp::row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, wv1);

    // get the differences between the new and the old columns
    pcps::agp::col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, wv2);

    // determine which column and row of the inverse matrix we need
    const int col = pcps::agp::find_index(nap, a_old_occ_sites[0], a_sites);
    const int row = pcps::agp::find_index(nap, b_old_occ_sites[0], b_sites);

    // special handling of the element at which the row and column cross
    wv1[row] = pcps::zero<SCALAR>();
    wv2[col] =   gamma[ a_new_occ_sites[0] * nas + b_new_occ_sites[0] ] - gamma[ a_old_occ_sites[0] * nas + b_old_occ_sites[0] ];

    // compute the action of the inverse matrix on the column difference
    pcps::xgemm('T', 'N', nap, 1, nap, pcps::unity<SCALAR>(), inv_mat, nap, wv2, nap,
                pcps::zero<SCALAR>(), wv3, nap);

    // compute the ratio of the new and old determinants
    retval =   ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + col, nap, wv1, 1) )
             * ( pcps::unity<SCALAR>() + pcps::xdot(nap, inv_mat + nap*row, 1, wv2, 1) )
             - inv_mat[ nap * row + col ] * pcps::xdot(nap, wv1, 1, wv3, 1);

    // check that the alpha and beta operators are arranged in one of the four allowed orders
    const int * const ds = diff_sites;
    assert( ( ds[0] % 2 == 0 && ds[1] % 2 == 0 && ds[2] % 2 == 1 && ds[3] % 2 == 1 ) ||    // RLRL = aabb
            ( ds[0] % 2 == 1 && ds[1] % 2 == 1 && ds[2] % 2 == 0 && ds[3] % 2 == 0 ) ||    // RLRL = bbaa
            ( ds[0] % 2 == 0 && ds[1] % 2 == 1 && ds[2] % 2 == 1 && ds[3] % 2 == 0 ) ||    // RLRL = abba
            ( ds[0] % 2 == 1 && ds[1] % 2 == 0 && ds[2] % 2 == 0 && ds[3] % 2 == 1 )    ); // RLRL = baab

    // if the operators are in the abba or baab ordering, apply a sign since we have assumed aabb ordering
    if ( ds[0] % 2 != ds[1] % 2 )
      retval = -retval;

  // unknown move type
  } else throw pcps::Exception( (boost::format("unknown move type in pcps::agp::det_ratio with nds = %i") % nds).str() );

  // return the determinant ratio
  return retval;

}

// explicitly instantiate the function
template double pcps::agp::det_ratio(const int nap,
                                     const int nas,
                                     const int nds,
                                     const int * const diff_sites,
                                     const int * const new_config,
                                     const int * const old_config,
                                     const double * const gamma,
                                     const double * const inv_mat,
                                     const int * const a_sites,
                                     const int * const b_sites,
                                     double * const wv1,
                                     double * const wv2,
                                     double * const wv3);
template std::complex<double> pcps::agp::det_ratio(const int nap,
                                                   const int nas,
                                                   const int nds,
                                                   const int * const diff_sites,
                                                   const int * const new_config,
                                                   const int * const old_config,
                                                   const std::complex<double> * const gamma,
                                                   const std::complex<double> * const inv_mat,
                                                   const int * const a_sites,
                                                   const int * const b_sites,
                                                   std::complex<double> * const wv1,
                                                   std::complex<double> * const wv2,
                                                   std::complex<double> * const wv3);

//-------------------------------------------------------------------------------
// pcps::agp::update -- Updates the inverse matrix and alpha/beta site lists for
//                      a new lattice configuration, which differs at the
//                      specified sites from the old configuration.
//-------------------------------------------------------------------------------

                                                            // exit status  array size    description
template <class SCALAR>
void pcps::agp::update(const int nap,                       // unchanged                  number of alpha particles
                       const int nas,                       // unchanged                  number of alpha lattice sites
                       const int nds,                       // unchanged                  the number of sites that differ between the configurations
                       const int * const diff_sites,        // unchanged      nds         the indices of the differing sites ( 0 <= index < 2*nas )
                       const int * const new_config,        // unchanged      2*nas       new lattice configuration
                       const int * const old_config,        // unchanged      2*nas       old lattice configuration
                       const SCALAR * const gamma,          // unchanged      nas*nas     the gamma matrix
                       SCALAR * const inv_mat,              // overwritten    nap*nap     the inverse of the agp matrix
                       int * const a_sites,                 // overwritten    nap         the occupied alpha sites ( numbered 0 to nas-1 )
                       int * const b_sites,                 // overwritten    nap         the occupied beta sites  ( numbered 0 to nas-1 )
                       SCALAR * const wv1,                  // overwritten    nap         work space
                       SCALAR * const wv2,                  // overwritten    nap         work space
                       SCALAR * const wv3)                  // overwritten    nap         work space
{

  // check that the number of sites affected by the move is supported
  assert(nds == 0 || nds == 2 || nds == 4);

  // there is nothing to do if no sites were affected
  if (nds == 0) return;

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int n_a_diff, n_b_diff;
  int a_new_occ_sites[2];
  int b_new_occ_sites[2];
  int a_old_occ_sites[2];
  int b_old_occ_sites[2];
  pcps::agp::get_diff_occs(nds, diff_sites, new_config, old_config, n_a_diff, n_b_diff,
                           a_new_occ_sites, b_new_occ_sites, a_old_occ_sites, b_old_occ_sites);

  // single alpha move
  if ( n_a_diff == 2 && n_b_diff == 0 ) {

    // apply a row update
    pcps::agp::row_update(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, a_sites, inv_mat, wv1, wv2, wv3);

  // single beta move
  } else if ( n_a_diff == 0 && n_b_diff == 2 ) {

    // apply a column update
    pcps::agp::col_update(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, b_sites, inv_mat, wv1, wv2, wv3);

  // double alpha move
  } else if ( n_a_diff == 4 && n_b_diff == 0 ) {

    // apply row updates
    pcps::agp::row_update(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, a_sites, inv_mat, wv1, wv2, wv3);
    pcps::agp::row_update(nap, nas, a_new_occ_sites[1], a_old_occ_sites[1], gamma, b_sites, a_sites, inv_mat, wv1, wv2, wv3);

  // double beta move
  } else if ( n_a_diff == 0 && n_b_diff == 4 ) {

    // apply column updates
    pcps::agp::col_update(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, b_sites, inv_mat, wv1, wv2, wv3);
    pcps::agp::col_update(nap, nas, b_new_occ_sites[1], b_old_occ_sites[1], gamma, a_sites, b_sites, inv_mat, wv1, wv2, wv3);

  // alpha beta move
  } else if ( n_a_diff == 2 && n_b_diff == 2 ) {

    // apply row update
    pcps::agp::row_update(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_sites, a_sites, inv_mat, wv1, wv2, wv3);

    // apply column update
    pcps::agp::col_update(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_sites, b_sites, inv_mat, wv1, wv2, wv3);

  // unknown move type
  } else throw pcps::Exception( (boost::format("unknown move type in pcps::agp::update with nds = %i") % nds).str() );

}

// explicitly instantiate the function
template void pcps::agp::update(const int nap,
                                const int nas,
                                const int nds,
                                const int * const diff_sites,
                                const int * const new_config,
                                const int * const old_config,
                                const double * const gamma,
                                double * const inv_mat,
                                int * const a_sites,
                                int * const b_sites,
                                double * const wv1,
                                double * const wv2,
                                double * const wv3);
template void pcps::agp::update(const int nap,
                                const int nas,
                                const int nds,
                                const int * const diff_sites,
                                const int * const new_config,
                                const int * const old_config,
                                const std::complex<double> * const gamma,
                                std::complex<double> * const inv_mat,
                                int * const a_sites,
                                int * const b_sites,
                                std::complex<double> * const wv1,
                                std::complex<double> * const wv2,
                                std::complex<double> * const wv3);

//-------------------------------------------------------------------------------
// pcps::agp::der_ratios -- Computes the ratios between the AGP determinant's
//                          derivatives and the AGP determinant's value for
//                          the requested variables.
//-------------------------------------------------------------------------------

                                                                   // exit status  array size    description
template <class SCALAR>
void pcps::agp::der_ratios(const pcps::Wavefunction<SCALAR> & wfn, // unchanged                  the wavefunction object
                           const int nap,                          // unchanged                  number of alpha particles
                           const int nas,                          // unchanged                  number of alpha lattice sites
                           const bool do_agp,                      // unchanged                  whether to compute derivatives for agp weights
                           const bool do_orb,                      // unchanged                  whether to compute derivatives for orbital coefficients
                           const bool do_par,                      // unchanged                  whether to compute derivatives for pairing matrix elements
                           const SCALAR * const gamma,             // unchanged      nas*nas     the gamma matrix
                           const SCALAR * const inv_mat,           // unchanged      nap*nap     the inverse of the agp matrix
                           const int * const a_sites,              // unchanged      nap         the occupied alpha sites ( numbered 0 to nas-1 )
                           const int * const b_sites,              // unchanged      nap         the occupied beta sites  ( numbered 0 to nas-1 )
                           SCALAR * const agp_der,                 // overwritten    nas         derivative ratios for agp weights
                           SCALAR * const orb_der,                 // overwritten    nas*nas     derivative ratios for orbital coefficients
                           SCALAR * const par_der)                 // overwritten    nas*nas     derivative ratios for pairing matrix elements
{

  // compute AGP wieght derivative ratios
  if (do_agp) {
    for (int i = 0; i < nas; i++) {
      SCALAR v = pcps::zero<SCALAR>();
      for (int row = 0; row < nap; row++) {
        v += pcps::unity<SCALAR>();
        for (int j = 0; j < nap; j++)
          v +=  inv_mat[j*nap+row]  *  (   wfn.restricted_orb_coeffs()[ i * nas + a_sites[row] ]
                                         * wfn.restricted_orb_coeffs()[ i * nas + b_sites[  j] ]
                                         - gamma[ a_sites[row] * nas + b_sites[j] ] );
        //     |                  |   |                                                          |
        //     --------------------   ------------------------------------------------------------
        //             /|\                                          /|\
        //              |                                            |
        //    matrix inverse element                difference between new and old row
      }
      agp_der[i] = v;
    }
  }

  // compute orbital coefficient derivative ratios
  if (do_orb) {

    pcps::xscal(nas * nas, pcps::zero<SCALAR>(), orb_der, 1);

    for (int i = 0; i < nas; i++) {

      // contribution from alpha particles
      for (int row = 0; row < nap; row++) {
        SCALAR v = pcps::unity<SCALAR>();
        for (int j = 0; j < nap; j++)
          v +=  inv_mat[j*nap+row]  *  (   wfn.agp_weights()[i] * wfn.restricted_orb_coeffs()[ i * nas + b_sites[  j] ]
                                         - gamma[ a_sites[row] * nas + b_sites[  j] ] );
        //     |                  |   |                                                                                |
        //     --------------------   ----------------------------------------------------------------------------------
        //                /|\                                              /|\
        //                 |                                                |
        //        matrix inverse element                     difference between new and old row
        orb_der[ i * nas + a_sites[row] ] += v;
      }

      // contribution from beta particles
      for (int col = 0; col < nap; col++) {
        SCALAR v = pcps::unity<SCALAR>();
        for (int j = 0; j < nap; j++)
          v +=  inv_mat[col*nap+j]  *  (   wfn.agp_weights()[i] * wfn.restricted_orb_coeffs()[ i * nas + a_sites[  j] ]
                                         - gamma[ a_sites[  j] * nas + b_sites[col] ] );
        //     |                  |   |                                                                                |
        //     --------------------   ----------------------------------------------------------------------------------
        //                /|\                                              /|\
        //                 |                                                |
        //        matrix inverse element                     difference between new and old column
        orb_der[ i * nas + b_sites[col] ] += v;
      }

    }

  }

  // compute pairing matrix element derivative ratios
  if (do_par) {
    pcps::xscal(nas * nas, pcps::zero<SCALAR>(), par_der, 1);
    for (int p = 0; p < nap; p++)
    for (int q = 0; q < nap; q++)
      par_der[ nas * a_sites[p] + b_sites[q] ] = inv_mat[ nap * q + p ];
  }

}

// explicitly instantiate the function
template void pcps::agp::der_ratios(const pcps::Wavefunction<double> & wfn,
                                    const int nap,
                                    const int nas,
                                    const bool do_agp,
                                    const bool do_orb,
                                    const bool do_par,
                                    const double * const gamma,
                                    const double * const inv_mat,
                                    const int * const a_sites,
                                    const int * const b_sites,
                                    double * const agp_der,
                                    double * const orb_der, 
                                    double * const par_der);
template void pcps::agp::der_ratios(const pcps::Wavefunction<std::complex<double> > & wfn,
                                    const int nap,
                                    const int nas,
                                    const bool do_agp,
                                    const bool do_orb,
                                    const bool do_par,
                                    const std::complex<double> * const gamma,
                                    const std::complex<double> * const inv_mat,
                                    const int * const a_sites,
                                    const int * const b_sites,
                                    std::complex<double> * const agp_der,
                                    std::complex<double> * const orb_der,
                                    std::complex<double> * const par_der);
